Loop updates for quantum Monte Carlo simulations in the canonical ensemble 
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We present a new non-local updating scheme for quantum Monte Carlo simulations, which conserves particle 
number and other symmetries. It allows exact symmetry projection and direct evaluation of the equal-time 
Green's function and other observables in the canonical ensemble. The method is applied to bosonic atoms in 
optical lattices, neutron pairs in atomic nuclei and electron pairs in ultrasmall superconducting grains. 
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Quantum Monte Carlo methods (QMC) allow an in princi- 
ple exact simulation of quantum many -body systems 1 1 ] . Over 
the past decade, cluster updates have increased the efficiency 
of lattice-QMC methods dramatically 1 2 , 3 ] . They proved par- 
ticularly useful near phase transitions, where the traditional 
algorithms suffered from a critical slowing down |4]. Based 
on an analogy with the Swendsen-Wang algorithm for classi- 
cal systems |5], the first loop algorithm for QMC constructed 
clusters in the form of loops, which then could be flipped 
to obtain new configurations \2]. An important development 
was the formulation in continuous imaginary time |6J], which 
eliminated discretization errors. The worm algorithm |7] es- 
tablished the link between the construction of the loops and 
the sampling in an extended configuration space, thereby al- 
lowing a direct evaluation of the one-body Green's function. 
The loop updates have also been implemented in the stochas- 
tic series expansion method (SSE) 1 8]. A further optimization 
of the loop construction process was developed in the form 
of directed loops \9]. These turned out to be a special case 
of a more general locally optimal strategy, which tries to op- 
timize the loop construction by using the optimal stochastic 
transition kernel at each local step of the process 
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forementioned loop algorithms all lead to fluctuations in par- 
ticle number (or in magnetization for spin systems). There- 
fore they sample the grand canonical ensemble. There are 
cases where particle-number conservation plays an important 
role and where one would like to sample directly the canon- 
ical ensemble: e.g. lattice systems at commensurate fill- 
ings | ll|] or finite systems such as ultrasmall superconduct- 
ing grains or atomic nuclei Results for the 
canonical ensemble can be obtained from loop algorithms by 
by using only those configurations in the sample which have 
the right particle number 1 16] or by explicitly throwing out the 
loop updates which change particle number 1 3]. However, this 
is impractical because a lot of effort is wasted on the discarded 
configurations or updates and because it requires a good esti- 
mate of the chemical potential. In this letter, we present a 
class of loop updates which explicitly conserve particle num- 
ber. The algorithm results in moves that are always accepted, 
which makes it easier to code and more efficient to run than 
other loop-update schemes. Furthermore, one can impose the 
conservation of other symmetries. 



Like most QMC methods, our algorithm starts from a 
decomposition of the imaginary time propagator, £/(p) = 
exp(— p//). Generally one can write the Hamiltonian as H — 
Hq — V, consisting of an easy part Ho and a residual interaction 
V (note the minus sign in front of V, in order to ease notations 
further on). For such a Hamiltonian, one can make a pertur- 
bative expansion in V using the following integral expression: 



f(P)=L JV(t l )V(t 2 )---V{t m )e-^dt l dt 2 ---dt m , (1) 



m=0 
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with V(t) = exp(—f//o)Vexp (f//o). Instead of sampling the 
representations of the operator U (p) directly, our method will 
perform a Markovian random walk in an extended config- 
uration space, related to the decomposition of the operator 
t/'(P,T) = e~ zH Ae~($~^ H , where A is the worm operator, to 
be defined later on. An alternative would be to insert a cre- 
ation and an anihilation operator at different imaginary times. 
This forms the basis of the continuous-time loop algorithm 
0] and the worm algorithm Q] . Here we will show that it is 
advantageous to work with a single worm operator, provided 
that it commutes with the residual interaction: AV — VA. If 
the operator A furthermore commutes with the generator of a 
symmetry of Ho and V, then one can restrict the configurations 
to specific symmetry representations, such that symmetry- 
projected results are obtained. In particular one can sample 
the canonical ensemble with a worm operator that conserves 
particle number. 

By taking the trace (restricted to the wanted particle num- 
ber and symmetry) and inserting complete sets of eigenstates 
of Hq between all operators in the integral representation of 
I/'(P,t), one assigns a weight W(m,i,t,z) to the configura- 
tions specified by an order m, a set of inserted eigenstates 
i'o , /i , . . . , i m , interaction times t \ , t% , . . . , t m , and the worm in- 
sertion time x. Let ii and ig denote the states to left and 
right of the worm operator. We will call the configurations 
for which il = ir diagonal configurations. One can choose 
the worm operator A such that its diagonal elements are con- 
stant, i.e. (i\A\i) = c for all basis states i. Then the sum of 
the weights of all the diagonal configurations is proportional 
to the particle-number projected trace of the operator t/(P), 
which is nothing else than the canonical partition function. 
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Hence, the sampling of the configurations proportional to their 
weight W(m,i,t,x) leads to a sampling of the canonical en- 
semble. 

To this end a Markov process is used for which the sta- 
tionary distribution reflects the contribution of the configu- 
rations to the extended partition sum Z' N ($) = Tr^v [C/'(P,t)]. 
Let us assume that a diagonal configuration is given. We are 
free to choose a new insertion time x for the worm opera- 
tor because the weight is independent of X. Then we per- 
form a number of Markov steps according to the following 
rules, until the move halts again in a diagonal configuration. 
We present the Markov rules in terms of a set of parameters 
Id, cd, ^Kdd 1 >£d, gD, a D and so, which are defined in detail in 
table|H 

• Evaluate the diagonal energy to the left (E£) and to the 
right (Er) of the worm operator, and evaluate the values 
of qi and qg. 

• Choose a direction D, either left (L) or right (R), propor- 
tional to the relative weights qi and qg. Let D' denote 
the opposite direction. 

• With probability cd, insert an interaction term V and 
choose a new intermediate state /, according to the dis- 
tribution P DD , (it) = (i D \A | U) (k\V | ijy ) I {i D \AV | i D < ) . 

• Move the worm in direction D over a step size Ax 
chosen according to a Poisson distribution: Ax 
exp(— EdAx). Assume periodic boundaries in imaginary 
time. 

• If the worm would pass an interaction term in the time 
step Ax, there are three options: (a) With probability 
so- remove the interaction and halt the worm. Then 
the Markov step ends here, (b) With probability ao- 
remove the interaction, adjust the parameters and con- 
tinue the worm move in the same direction, (c) With 
probability 1 — ap — so, or if the interaction can not be 
removed: let the worm pass the interaction, but choose 
a new intermediate state according to Poiyi'd- Adjust 
the parameters and continue the worm move in the same 
direction. 

• When the worm has reached the end of the time step 
Ax, then with probability 1 — go the worm is halted and 
the Markov step ends here. With probability go, insert 
an interaction and choose a new intermediate state ac- 
cording to P DD i(ij). Adjust the parameters, draw a new 
time step Ax and continue the worm move in the same 
direction. 

There is considerable freedom in the definition of the sam- 
pling parameters. Guided by the principle of locally opti- 
mal moves 11(1 Il7ll . we propose the definitions as given in 
table H] These rules assure that the Markov chain satisfies 
the detailed balance condition for the weight Ru{W(m,i,t,z), 
with Rir = qg + qL- One could use the Metropolis-Hastings 



algorithm 1 18l ll9ll in order to sample according to the weights 
W(m, i,t,x). Because the factors Rlr fluctuate only mildly in 
practice, one can accept all moves and take the extra weight- 
ing factor Rlr into account when evaluating observables. This 
speeds up the algorithm and reduces the complexity of the 
code. Note that during a Markov step the worm keeps moving 
in the same direction, even after passing, inserting or remov- 
ing interactions. This is possible while maintaining detailed 
balance because A and V commute. Otherwise one would 
have to consider so-called bounces: the worm could return 
on the path were it came from and undo its last changes. This 
is one of the reasons why our new algorithm is more efficient 
than the worm and loop algorithms of Refs. J3, 00 El- 
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TABLE I: Definitions for the sampling parameters, assuming that 
Ed < Ed (interchange D and D 1 otherwise). The global parameter 
<|) should be taken small (such that < 9{dd for all diagonal config- 
urations) but not zero, to make sure that the worm halts in diagonal 
configurations that are sufficiently decorrelated. 

While each Markov step is based on local changes, the 
chain of steps between two diagonal configurations corre- 
sponds to a global loop update: if one follows the creation and 
annihilation part of the worm operator as they move through 
the world-line representation of the configurations, one sees 
that they describe a loop which closes again when the con- 
figuration becomes diagonal. By keeping track of the worm 
operator at the intermediate steps, one can collect statistics 
for the expectation values of non-diagonal operators, similar 
to the way one evaluates the one-body Green's function in the 
worm algorithm |7]. The advantage of our method is that the 
operators are always evaluated at equal imaginary times, lead- 
ing to much better statistics for non-diagonal operators. 

As an illustration we have applied our method to the 
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one-dimensional Bose-Hubbard Hamiltonian, given by H = 
— *X;&Jfc;+i +t/Li n i( n i — 1) /2, with rii — b]bj. It models 
the low-energy degrees of freedom of various physical sys- 
tems: cold bosonic atoms in an optical lattice |20], 4 He atoms 
on graphite |21], and superconducting islands or grains con- 
nected by Josephson junctions 1 22] . Various many -body tech- 
niques have been deployed to study its phase diagram. With- 
out trying to be complete, we mention here algebraic 12 311 . 
mean-field [11], perturbative I24L renormalization- group ll25ll 
and QMC approaches lEfill . For our method we take V = 
Y^ibjbi+i an d A — Y,i.jb\bj. At commensurate fillings, this 
model can undergo a quantum phase transition from a super- 
fluid to a Mott isolator. With grand-canonical methods it is 
difficult to simulate this transition because one has to deter- 
mine the chemical potential such that the exact density is ob- 
tained. Fig. ^ shows the superfluid and condensed fraction 
for a uniform one-dimensional system of 128 sites at nearly 
zero temperature (p = T~ l = 128f~') at a density of exactly 
one particle per site. The superfluid fraction was derived from 
the winding numbers [27], while the condensed fraction was 
obtained from the one-body equal-time Green's function. 
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FIG. 1: Superfluid (p. s ) and condensed fraction (p c ) for the one- 
dimensional Bose-Hubbard model on a uniform lattice of 128 sites, 
at an inverse temperature P = 128f~' . 



Fermionic models can be handled too, provided that there 
is a symmetry which guarantees that all matrix elements are 
positive, such that the sign problem is absent. This is the case 
for the pairing Hamiltonian: 



H 



E 

a jm 



e j a ajm a aj™ 



E 

a jma'j'm' 



J ajm u a/m"<x' /";«'"</.' i' in' ■ 



(2) 

where m creates a particle in a state with quantum numbers 

a,j,m, and a^j m a particle in the time-conjugated state. The 
pairing model has been used extensively in nuclear physics 
to model the pair correlations between nucleons and to ex- 
plain the particular odd-even effects found in nuclear spec- 
troscopy b(Sl I29I1 - Canonical QMC methods for the pair- 
ing Hamiltonian have been presented before: Cerf's world- 



line QMC method 1 30] does not have a sign problem, but it 
does not sample the full phase space of broken pairs l3lll . 
which makes it impractical for finite temperature calculations. 
The shell model Monte Carlo method overcomes this prob- 
lem II4II15LI32II . but has a sign problem for odd particle num- 
bers. For a constant pairing strength G, the eigenstates can be 
calculated exactly 13311341 . Our method supplements these al- 
gebraic solutions with finite temperature results. It can also be 
applied to level-dependent pairing interactions, for which no 
algebraic solution is available. Here, we take V to be the non- 
diagonal part of the pairing interaction and we take A equal 
to V plus a constant diagonal term. A pair-breaking term has 
to be added to A and V in order to ensure ergodicity fnll . 
Note that the worm operator conserves angular momentum. 
Therefore one can restrict the intermediate states to a specific 
value of the quantum numbers J and J z . Fig. [2] shows how 
the .^-projected finite-temperature results converge to the ex- 
act eigenvalues at low temperature, for a model that describes 
neutron pair correlations in 56 Fe 11511 . 
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FIG. 2: /--projected internal energies as a function of temperature 
for 10 neutrons in the pf + sdg shell, for a constant pairing interac- 
tion with G = 0.286 and single-particle energies as in Ref. 11511 . On 
the left hand side the exact eigenvalues are shown 1 34]. 



The pairing model where all levels are double degener- 
ate ( j =1/2 for all a) and equally spaced (with level spac- 
ing d), is called the picket fence model. It applies to ultra- 
small metallic grains, provided that a canonical approach is 
used fl2l IT3I l35ll . One has found that grand-canonical ap- 
proaches lead to an abrupt but unphysical suppression of the 
superconductive correlations for level spacings larger than the 
bulk gap A. The extension of the exact ground-state cal- 
culations to finite temperatures was cited as an open prob- 
lem in Ref. 1 13], and has only been treated at the mean-field 
level llirjl (except for the smallest systems). Our method pro- 
vides exact finite-temperature results both for odd and even 
particle numbers. Therefore we can study the odd-even asym- 
metry, which is a key indicator of superconductive correla- 
tions. An illustrative quantity is the canonical pairing gap 
A can , as defined by Eq.(92) of Ref. 11 311 . It is a measure of 
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the change in energy due to pairing correlations. Its deviation 
from the BCS bulk gap indicates the difference between the 
canonical and grand-canonical ensemble. In Fig. [3] one sees 
the canonical pairing gap converging to the bulk gap with in- 
creasing system size. There is a clear odd-even effect at lower 
temperatures that disappears at higher temperatures, and this 
transition temperature decreases with increasing system size. 




T/A 

FIG. 3: Canonical pairing gap A can as a function of temperature, in 
units of the bulk gap A, for the picket fence model at half filling (N 
particles), for various system sizes (44, 86, 174 and 300 levels) and 
with pairing interaction strength G = 0.224d (as in Ref. |35]). 



These applications show the versatileness of the canoni- 
cal loop updates. Our method can be used to sample con- 
figurations with specific symmetries and in particular to sam- 
ple the canonical ensemble. It leads to a very efficient sam- 
pling scheme with all moves accepted and without 'bounces' 
or critical slowing down. It can be applied to bosons and 
fermion pairs in discrete model spaces and to spin systems 
at fixed magnetization. Off-diagonal observables such as the 
equal-time one-body Green's function can be evaluated with 
high efficiency. This opens new perspectives for the study of 
quantum many -body systems where particle number and other 
symmetries play an important role, such as the nuclear shell 
model and ultrasmall superconducting grains. 
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